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Abstract 


Phase separation mechanisms can produce a variety of complicated and intricate microstructures, which 
often can be difficult to characterize in a quantitative way. In recent years, a number of novel topological 
metrics for microstructures have been proposed, which measure essential connectivity information and are 
based on techniques from algebraic topology. Such metrics are inherently computable using computational 
homology, provided the microstructures are discretized using a thresholding process. However, while in 
many cases the thresholding is straightforward, noise and measurement errors can lead to misleading 
metric values. In such situations, persistence landscapes have been proposed as a natural topology metric. 
Common to all of these approaches is the enormous data reduction, which passes from complicated 
patterns to discrete information. It is therefore natural to wonder what type of information is actually 
retained by the topology. In the present paper, we demonstrate that averaged persistence landscapes 
can be used to recover central system information in the Cahn-Hilliard theory of phase separation. More 
precisely, we show that topological information of evolving microstructures alone suffices to accurately 
detect both concentration information and the actual decomposition stage of a data snapshot. Considering 
that persistent homology only measures discrete connectivity information, regardless of the size of the 
topological features, these results indicate that the system parameters in a phase separation process affect 
the topology considerably more than anticipated. We believe that the methods discussed in this paper 
could provide a valuable tool for relating experimental data to model simulations. 

1 Introduction 

Complicated patterns which evolve with time occur in a variety of applied contexts, and quantifying or even 
just distinguishing such patterns can pose serious challenges. Over the last decade, computational topology 
has emerged as a tool which on the one hand allows for significant data reduction, while at the same 
time focusing on and retaining essential connectivity information of the studied patterns. One particularly 
interesting application area is materials science, where complex evolving patterns are frequently created 
through complicated phase separation processes. Computational topology encompasses a wide variety of 
possible tools CH], and most of them in one way or another are based on homology theory. The deeper 
reason for this can be found in the inherent computability of homology, and in recent years powerful algorithms 
have been devised which enable fast homology computations for very large data sets, see for example nans], 
as well as the references therein. 

Homology has been used in a number of materials science contexts. Earlier studies have made use of the 
easily computable Euler characteristic, see for example the references in the recent survey [3l|. However, it 
has been pointed out that in some situations the information retained by the Euler characteristic is not enough 
to distinguish certain important pattern features. In contrast, the Betti numbers, which are associated with 
the homology groups and will be described below, provide a finer metric. They were used for example in CZI 
to relate the pattern complexity evolution as described by averaged Betti number evolution curves to the 
amount of stochasticity inherent in a phase field model for binary phase separation in metal alloys. This 
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study was the first to use homology information in the context of model validation. Based on available 
experimental data it was shown that if the noise in the system is too low, the observed Betti number 
evolution curves are qualitatively different from the experimental ones. In addition, it was demonstrated 
in nzi that while Betti numbers can be used to separate bulk from boundary behavior, the averaged Euler 
characteristic can only describe the boundary effects. Similar materials science studies in the context of 
polycrystals can be found in (2S1135]. While the first of these papers uses homology to study the connectivity 
properties of grain-boundary networks in planar sections of polycrystals, the second paper employs Betti 
numbers as a means to describe the thermal-elastic response of calcite-based polycrystalline materials such 
as marble. More precisely, in m homological techniques are used to characterize not only the elastic energy 
density and maximum principal stress response fields in a polycrystal, but also the respective grain-boundary 
misorientation distributions which generated these response fields. It was shown that this topological analysis 
can quantitatively distinguish between different types of grain-boundary misorientations, and relate them to 
differences in the resulting response fields. 

In all of the applications described so far, the numerical or experimental data is given in the form of a 
field, or in other words, a real-valued function t/ : Q —)■ M defined on some domain Q c The associated 
patterns are subsets of the domain Q, and usually created through a thresholding process. For example, 
after selecting a suitable threshold 6, one can consider the sub- and super-level sets 

= {x e Q : ± (t/(x) - 0) > 0} , 

which in some contexts are called nodal domains of u. As subsets of Q c the sets have well- 
defined singular homology groups, and if the function u is sufficiently regular, these groups can be computed 
precisely; see for example dKia for the two-dimensional case d = 2. However, if the function u is not 
smooth, or if the thresholding process involves a field created from experimental or noisy data, then the above 
thresholding process may not capture the correct topology of the actual underlying pattern. Moreover, in 
certain applications the thresholding approach itself might not be appropriate, for example if there is no 
obvious or physically relevant choice of threshold. 

An extension of the concept of homology to situations involving noise or the lack of a clear thresholding 
process has been proposed some fifteen years ago and is called persistent homology [18]. While this extension 
is described in more detail in the next section, persistent homology is a dimension reduction technique which 
provides a topological description of the evolution of sublevel sets of a mapping u as a function of the 
thresholding level 6. It gives rise to intervals [0i,02] over which certain topological features persist, and 
the length 02 — 0i can in some sense be viewed as a measure of importance of the specific feature. Thus, 
topological features with small interval lengths are usually considered as “noise" or "unimportant,” while 
features with long intervals are deemed “significant." Needless to say, the precise meaning of these notions 
will change from application to application. The concept of persistent homology has already been used in 
a number of materials science contexts, such as for example in the analysis of granular media |2Ql |2I], in 
the study of protein compressibility [T6|, as well as in the classification of amorphous structures m and 
glass I2HI- Efficient algorithms for computing persistent homology are described in |TH11^ . 

Despite the success of the above uses of computational topology in applications, one immediate question is 
the extent of the resulting dimension reduction. Through homology, patterns or microstructures are basically 
reduced to a finite set of integers, and it is therefore natural to wonder what information is still encoded in 
this reduced measurement, beyond the obvious connectivity information. It was pointed out in m that even 
small Betti number counts, which taken in isolation do not provide much in terms of pattern differentiation, 
can provide significant information when viewed in a stochastic, i.e., averaged setting. More precisely, it was 
shown in m that during the phase separation process called nucleation, averaged droplet counts on small 
domains give extremely precise projections for the observed averaged droplet counts on large domains, where 
the droplet count is Just the zero-dimensional Betti number. 

In the present paper, we demonstrate that when viewed in a stochastic and time evolving framework, 
topological information encodes considerably more than anticipated. This will be accomplished in the setting 
of instantaneous phase separation in binary metal alloys, as modelled by the Cahn-Hilliard theory of spinodal 
decomposition 019]. This phase separation phenomenon is initiated immediately after a high-temperature 
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Figure 1: Phase separation through spinodal decomposition in the Cahn-Flilliard-Cook model @ with white 
noise forcing. The images show snapshots of solutions which originate at the homogeneous state u = fj,, 
and for parameter values e = 0.005 and a = 0.001. From top left to bottom right the images are for 
mass jj, = 0.07 • k, where /c = 0,..., 7. 


melt of two uniformly mixed metal components is quenched, i.e., rapidly cooled. Depending on the concen¬ 
trations of the involved components, they will quickly separate to form complicated microstructures which 
contain some apparent element of randomness. Some of the resulting patterns in two space dimensions 
are shown in Figure These patterns evolve with time, and after the initial phase separation process a 
coarsening stage sets in, during which the characteristic length scale of the patterns increases. 

The first mathematical model for spinodal decomposition was introduced by Cahn and Hilliard izmo], 
who proposed a nonlinear evolution equation for the relative concentration difference u = pa — Pb, where pA 
and Pb denote the relative concentrations of the two components, i.e., Pa + Pb = 1- Their model is based 
on the Ginzburg-Landau free energy given by 

Ee{u) = l^(^^\Vuf + F{u)^ dx, (1) 

where Q c is a bounded domain, and the positive parameter e models interaction distance. The bulk 
free energy F is a double well potential, which for the purposes of this paper is taken as 

F(u) = 1)" . (2) 

Taking the variational derivative 5E^/5u of the Ginzburg-Landau free energy © with respect to the concen¬ 
tration variable u, one then obtains first the chemical potential 1 / 1 / = —e^A^y-h F'{u), and then the associated 
Cahn-Hilliard equation dujdt = Aw, i.e., the fourth-order partial differential equation 

^ = -A{e^Au-F'{u)) , (3) 

subject to homogeneous Neumann boundary conditions for both w and u. Due to these boundary conditions, 
any mass flux through the boundary is prohibited, and therefore mass is conserved. We generally consider 
initial conditions for Q which are small-amplitude random perturbations of a spatially homogeneous state, 
i.e., we assume that u{0, x) th p, for all x e Q, as well as u(0, x) dx/\Q\ = p, where we use the standard 
abbreviation |Q| = JqI dx. One can easily see that such initial conditions lead to instantaneous phase 
separation in the Cahn-Hilliard equation as long as the conserved total mass p satisfies F"{p) < 0. In the 
context of ©, this condition is equivalent to \p\ < \l\fz. 
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While the deterministic Cahn-Hilliard model ([^ provides a basic qualitative model for spinodal decomposi¬ 
tion, it cannot produce microstructures whose evolution during the initial phase separation stage agrees with 
experiments m As it turns out, this deficiency is due to the fact that the original model @ completely 
ignores thermal fluctuations. To remedy this, Cook [12] extended the model by adding a random fluctuation 
term £, i.e., he considered the stochastic Cahn-Hilliard-Cook model 

f) n 

— = (4) 

where the expected value and the correlation of the noise process satisfy both 

(^(t,x)) = 0 and Xi)^(t 2 , X 2 )) = 5(fi - t 2 )c 7 (xi - X 2 ) . 

Here a > 0 is a measure for the intensity of the fluctuation and q describes the spatial correlation of the 
noise. Thus, the noise is always uncorrelated in time, and for the special case q = 5 vje obtain space-time 
white noise. In general, however, the noise process will exhibit spatial correlations. For more details we 
refer the reader to the discussions in (HI [22] • As mentioned before, while both the deterministic and the 
stochastic model generate microstructures which are qualitatively similar to the ones observed during spinodal 
decomposition [H], only the stochastic model can in principle lead to a more quantitative agreement m 
Recent mathematical results for the models © and @ have identified the observed microstructures as 
certain random superpositions of eigenfunctions of the Laplacian, and were able to explain the dynamics of 
the decomposition process in more detail [2l l4l (3011^ l33] . 

Based on the insight obtained in Cz], our studies in the present paper concentrate exclusively on the 
stochastic Cahn-Hilliard-Cook model ©• For this model, we show that persistent homology averages suffice 
to accurately deduce the total mass q,, and therefore the total alloy component concentrations, as well as the 
actual stage of the decomposition process. Both of these observations are somewhat surprising, since there is 
no obvious link between the connectivity information retained by homology and either of these quantities. Our 
results demonstrate, however, that during the phase separation process, the evolution of averaged persistent 
homology information in some sense is a very detailed encoding of these system parameters. Expressed 
differently, while each mass value q uniquely determines an averaged persistent homology evolution for 
randomly selected solutions originating close to the homogeneous state q, the topology encoded in these 
averaged evolutions is sufficiently different to allow for the solution of an inverse problem. A similar statement 
can be made for the stage of the decomposition. Moreover, we wold like to point out that unlike the Betti 
numbers, persistent homology information cannot easily be averaged directly. We therefore make use of the 
recent concept of persistence landscapes EIE], which will be described in more detail below. In our studies, 
we are also using the idea of topological process, which will be introduced later in this paper. It allows us to 
deal with a processes that produces a sequence of topological spaces. 

The remainder of this paper is organized as follows. In Section]^ we survey the topological background 
for this paper, to keep our presentation self-contained and accessible for applied scientists. After introducing 
cubical homology in Section |2.1| we discuss the concept of persistent homology in Section |2.2[ Since 
for our study the averaging of persistence information is crucial, we then turn our attention to persistence 
landscapes in Section After these preparations, the topological analysis of microstructures created by the 
Cahn-Hilliard-Cook model © is the subject of Section We begin by presenting the basic methodology in 
Section |TI] turn our attention to the topological determination of the total mass q in Section and show 
in Section [T^ that even the time of a solution snapshot can be recovered accurately from the persistence 
information during the spinodal decomposition and the initial coarsening regime. Finally, Section [^contains 
a brief summary and draws some conclusions. 

2 Topological Background 

In order to keep this paper as self-contained as possible, this section provides background information on 
homology and persistence. More precisely, we recall the basic definition of cubical homology, introduce 
persistent homology, and finally describe the concept of persistence landscapes. 
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Figure 2: Illustration of two fundamental concepts from cubical homology, (a) Boundary Operator: The 
left image shows the two cubes Ci = [1,1] x [1,2] and C 2 = [2,3] x [2,3], whose dimensions are one and 
two, respectively. Their boundaries can easily be computed as dCi = [1,1] x [1,1] + [1,1] x [2,2], as well 
as dC 2 = [2, 2] X [2, 3] + [3, 3] x [2,3] + [2, 3] x [2, 2] + [2,3] x [3, 3]. Note that in both cases, the formal 
notion of boundary is consistent with the intuitive boundary of that cube, (b) Cycles: Consider the cubical 
complex shown in gray in the right image, together with all contained edges and vertices. Then both the 
collection of edges shown in red, and the collection shown in green form one-dimensional chains which are 
cycles, i.e., their boundaries vanish. This is due to the fact that every cell of dimension zero in the support 
of each chain occurs in the boundary of exactly two cells, and therefore cancels out in the chain's boundary. 


2.1 Cubical Homology 

Homology theory associates discrete objects with topological spaces in such a way that continuous defor¬ 
mations of the space do not change the objects. In other words, homology is a topological invariant of the 
underlying topological space. We will see below that the discrete objects are Abelian groups, and depending 
on the form in which the topological space is given, a number of different homology theories have been 
developed over the years to compute them. If a given space can be treated within several different homology 
theories, then the homology groups assigned by these different homology theories will be isomorphic. For the 
purposes of this paper, we restrict our attention to one specific homology theory, namely cubical homology. 
For this, we assume that our topological space is given as a specific union of cubical sets of possibly varying 
dimensions, leading to the concept of a cubical complex. 

To be more precise, define an elementary interval as a compact real interval of the type [i,i+ 1] or 
where I denotes an integer. While the first interval type is referred to as non-degenerate, the second one is 
called a degenerate interval. An elementary cube C is a product of elementary intervals /i x ... x /„, and its 
dimension equals the number of non-degenerate intervals 1^ in this product. Furthermore, the total number n 
of intervals in the product is called the embedding dimension of the cube C. For example, Figure|^aJ depicts 
an elementary cube of dimension one in black, and a two-dimensional elementary cube in blue, and both cubes 
have embedding dimension two. 

While elementary cubes form the building blocks of the cubical complexes introduced below, we first need 
to introduce a few algebraic concepts based on them. A chain c is a formal sum 

c = ctiCi -F ... -F ocinCin , where Q. Cm are elementary cubes and cti. € F are scalars. (5) 

In this formula, F denotes an arbitrary field. However, for the purposes of this paper, we always consider the 
field F = Z 2 which consists of the two elements 0 and 1, as well as addition and multiplication modulo two. 
Thus, for the chain c in Q a cube Q is either present or absent, depending on whether = 1 or = 0, 
respectively. Furthermore, the support of the chain c is the collection of all elementary cubes Q which 
have nonzero coefficient ai^ in c. Two one-dimensional chains consisting of four one-dimensional elementary 
cubes each are shown in green and red in Figure 

Next we have to introduce the algebraic boundary of a chain. For an elementary interval / = [n,m], 
its boundary is the chain supported by the two degenerate intervals [m, m] and [n, n], and this immediately 
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implies both 


as well as 


5[n, /? + 1] = [/? + 1, /? + 1] + [n, n] 


d[n, n] = 0 , 


depending on whether m = n + l or m = n, respectively. Note that if the interval / is degenerate, we have 
used the fact that d[n, n] = [n, n] + [n, n] = (1 +!)[/?, n] = 0 in the field F = Z 2 . Combining the definition of 
the boundary for elementary intervals with the definition of elementary cubes, one can then lift the boundary 
definition to the latter using linearity. More precisely, let C = /i x ... x /„ denote an elementary cube. Then 
we define 

n 

dc = d{li X ... X In) = ^2 /l X ■ ■ ■ X 9// X ■ ■ ■ X /n . 

i=l 

where /i x ... x d[m, m + 1] x ... x F = li x ... x [m + 1, m + 1] x ... x !„ + li x ... x [m, m] x ... x \r\ 
the nondegenerate case, and li x ... x d[m, m] x ... x /„ = 0 in the degenerate one. Finally, the boundary 
operator extends linearly from elementary cubes to chains, i.e., for a chain c = cciCi +... Tcc^C^ one defines 
its boundary via dc = aidCi + ... + a^dCi. These concepts are illustrated further in Figure 

The boundary operator formalizes the intuitive notion of the boundary of a /c-dimensional elementary cube. 
For k = 0, such a cube is a point, and its boundary is zero, while for /c = 1 the cube is an interval with the 
boundary consisting of its endpoints. Note, however, that this notion does not depend on the embedding 
dimension of the cube, in contrast to the topological boundary. In other words, the boundary dC is determined 
from the intrinsic dimension of the cube. This latter fact also lies at the heart of the fundamental property 
of the boundary operator, which states that 


ddc = 0 for every chain c . 


One can easily see why this statement is true for elementary cubes. If C is an elementary cube of dimension k, 
then its boundary is a union of elementary cubes of dimension k — 1. This implies that the chain ddC consists 
of cubes of dimension k — 2. Every one of these k — 2-dimensional cubes is contained in the boundary of 
exactly two k — 1-dimensional cubes in dC, i.e., the chain ddC contains every elementary cube exactly 
twice. Since we are considering the field F = Z 2 , these cubes cancel out and yield a zero boundary. This is 
illustrated in Figure 

As we mentioned earlier, elementary cubes are the fundamental building blocks for the topological spaces 
which can be treated using cubical homology. Loosely speaking, we will consider subsets of E" which can 
be obtained as unions of elementary cubes of embedding dimension n. More precisely, we say that a finite 
collection /C of elementary cubes with embedding dimension /? is a cubical complex, if for every cube C € /C 
all elementary cubes in its boundary dC are also contained in K.. For a given cubical complex, we can then 
talk about its associated chains. For this, let 0 < k < n be fixed. Then chains which consist exclusively 
of elementary cubes in 1C with dimension k are called k-chains of 1C, and the set of all k-chains is denoted 
by Ck{IC). One can easily see that Q(/C) is an additive group, where again we use coefficients in F = Z 2 , and 
we will therefore refer to Ck{IC) as the k'^^ chain group. Furthermore, if the boundary operator d is applied 
to a k-chain c, then dc is clearly a k — 1-chain. In fact, one can show that the map d : Q(/C) ^ Q_i(/C) 
is a group homomorphism. Since we consider coefficients in the field F = Z 2 , even more is true: The chain 
groups Q(/C) are vector spaces over Z 2 , and the maps d : Q(/C) ^ Q_i(/C) are all linear. 

We are now finally in a position to describe cubical homology. For this, let K. denote a cubical complex. 
Intuitively speaking, homology is used to count the number of k-dimensional holes in the cubical complex, 
for every 0 < k < n. As a first step towards this goal, we call a k-chain c G Ck{IC) a k-cycle, if its boundary 
vanishes, i.e., if we have dc = 0. The set of all k-cycles is denoted by Zk{IC), and one can easily see that it 
is a subgroup of the k^^ chain group. For example, in Figuretwo 1-cycles are shown in red and green. 
It is clear from this image, that cycles alone cannot be used to detect holes in the cubical complex 1C. While 
the red 1-cycle in the figure encloses a hole, the green one does not. 

In order to distinguish between cycles which enclose holes and cycles which do not, we need to introduce 
one final concept. For this, assume as before that /C is a cubical complex and that 0 < k < n. Then a 
k-cycle is called a k-boundary, if it is equal to the boundary of some k -F 1-chain. We denote the set of 
all k-boundaries by Bk{IC), and since this set forms a subgroup of the k^'^ chain group, it is called the k^'^ 
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boundary group. In fact, due to the fundamental property 99 = 0 of the boundary operator, one has the 
inclusion Bj^{lC) c Zj^{lC), i.e., the /c-boundaries form a subgroup of the /c-cycles. We now return to the two 
1-cycles shown in Figure Notice that the green cycle is the boundary of the 2-chain which consists of 
all elementary cubes in the interior of the green loop. In other words, the green cycle is bounding a region 
of 2-dimensional cubes in the cubical complex, which completely fill its interior. On the other hand, one can 
show that the red 1-cycle cannot be represented as the boundary of a 2-chain, and this cycle does enclose a 
hole in the cubical set. 

The last example captures the essence of homology. In order to detect holes in a cubical complex, one 
needs to consider /c-cycles, as they are the /c-chains without boundary. On the other hand, one needs to 
make sure that such a /c-cycle actually encloses a hole, and for this one has to make sure that the cycle is 
not the boundary of a /c -F 1-chain in 1C. Otherwise, this k + 1-chain would fill up the holes created by the 
cycle. Notice, however, that just counting the /c-cycles which are not /c-boundaries does not give the correct 
number of holes. In fact, in Figure one can easily find 1-cycles which are different from the red one, 
but which still enclose the same hole. While not quite as intuitive, one can show that two cycles enclose the 
same holes, if their difference happens to be a /c-boundary. 

From a mathematical point of view, counting the number of holes seems to be related to counting equiv¬ 
alence classes. Call two /c-cycles homologous, if their difference is contained in the boundary group 6/<(/C). 
This defines an equivalence relation on the set of all /c-cycles. Furthermore, a cycle is is homologous to zero 
if and only if it is a /c-boundary, i.e., if it does not enclose a hole. Based on this discussion, it should not 
come as a surprise that the quotient groups 

Hk{}C) = Zk{}C)/Bk{IC) for all 0<k<n 

capture information about /c-dimensional holes in the cubical complex 1C. These groups are called homology 
groups, and since we consider coefficients in the field F = Z 2 , they are in fact quotient vector spaces. 

The homology groups defined above are abstract algebraic objects which can be assigned to a cubical 
complex. But what exactly is their relation to the number of /c-dimensional holes in fC. To see this, notice 
that if we have two /c-cycles ci and C 2 which enclose to different holes in the cubical complex, then the 
equivalence class in Hk{}C) determined by their sum ci + C 2 encloses both holes. In other words, two 
different holes in K, give rise to three different equivalence classes in the k^'^ homology group. In fact, what 
we need are the independent classes that can be found in Hk{IC), i.e., we need to determine the size of a 
basis in the quotient vector space Hk{)C), which is the same as the rank of the quotient group. If we now 
define 

/3/((/C) = rank Hk{IC) , called the k^^ Betti number of fC , 

then the number of /c-dimensional holes in the cubical complex is given by Pk{lC). Furthermore, the k^^ 
homology group is isomorphic to Intuitively, the zero-dimensional Betti number counts the number 

of connected components of the cubical complex, while the one-dimensional Betti number counts the number 
of tunnels or one-dimensional holes. Similarly, the two-dimensional Betti number counts the number of voids 
or cavities in the complex, and this classification of holes continues for higher dimensions. For example, the 
cubical complex in Figure is connected and has one hole, but as a flat object no cavities. This implies 
both /3o(/C) = 1 and /3i(/C) = 1, while /32(/C) = 0. 

The cubical homology theory introduced above provides a very rough quantitative measurement for the 
topological complexity of a complex fC and its underlying topological space. To make this more precise, 
let |/C| C M" denote the union of all elementary cubes in K,. Then one can show that the information provided 
by the homology groups is invariant under continuous transformations of |/C|. Specifically, if two cubical 
complexes ICi and IC 2 are such that their underlying topological spaces |/Ci| and I/C 2 I are homeomorphic, 
then the homology groups Hi^{JCi) and /-//((/C 2 ) are isomorphic. In fact, this is still true even if the spaces |/Ci| 
and I/C 2 I are only homotopy equivalent. In other words, homology captures connectivity and hole information 
that is invariant under continuous deformations. 
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Figure 3: Cubical complexes and filtration assigned to a gray-scale image. The image is shown in (a) and 
consists of sixteen pixels. From this image, one can extract cubical subcomplexes by considering only pixels 
with gray-scale intensity less than a given threshold. Depending on the threshold, this leads to one of the 
five cubical complexes shown in (b). Alternatively, one can assign one filtration of cubical complexes to the 
gray-scale image, as shown in (c). This filtration consists of five nested cubical complexes. 


2.2 Persistent Homology 

The homology groups introduced in the last section are topological invariants which can be assigned to a 
fixed topological space |/C|. In applications, finding the “correct" space to study is often a straightforward 
task, but sometimes the situation is less clear. Consider for example a gray-scale image, and suppose we are 
trying to use homology to detect holes in the image. One simple instance of this is shown in Figure|^a^. In 
this 4 X 4-pixel image, the lighter-colored pixels seem to enclose a hole given by the solid black pixel. If we 
are trying to use homology to detect this hole, we need to create a cubical complex which consists of the 
"correct pixels" of the image. This can be accomplished for example by choosing a gray-scale threshold, and 
then defining K. to consist of all pixels with intensity less than that threshold, together with all their edges 
and vertices. Depending on the threshold, one obtains the cubical complexes shown in Figureand the 
third and fourth of these complexes do correctly identify the hole visible in the image. Note, however, that 
the remaining three complexes fail to detect the hole. 

But what about a situation in which the image contains two holes, whose interiors occur at different 
gray-scale levels? Depending on the intensities of the respective surrounding pixels, it is conceivable that one 
cannot find a threshold value which detects both holes at the same time. In situations like this, finding one 
threshold level for the whole image clearly is the wrong goal. Rather, it would be nice to have an automatic 
way which allows for the detection of holes at different levels — and using homology in its original form 
alone makes this a difficult task. 

To address problems of the above type, the concept of persistent homology has been proposed, and it 
will be described in more detail in the following. For a literature review on the subject, we refer the reader 
to [T5] and the references therein. The fundamental idea behind persistent homology can easily be motivated 
by the complexes shown in Figure]^ Rather than singling out a specific cubical complex in (b), persistent 
homology deals with all of these complexes using the concept of a filtration. In the situation of the figure, the 
underlying filtration is the increasing (with respect to set inclusion) sequence of the five cubical complexes 
shown in (c), indexed by the threshold levels indicated by (1) through (5). In other words, a filtration encodes 
how the full image is assembled by adding pixels according to their gray-scale level, where pixels of the same 
level will be added simultaneously. Even though we only mention the pixels, which are two-dimensional cubes, 
we implicitly assume that at every stage of the filtration we do have cubical complexes as defined in the last 
section. This means that the edges and vertices of the pixels have to be added at the same time as the pixel. 
This is indicated in Figure|^cj by highlighting the one- and zero-dimensional cubes. This implies for example 
that at stage (1) the shown cubical complex has four vertices, four edges, and one square, while at stage (2) 
we have eight vertices, ten edges, and three squares. The goal of persistent homology is to detect and keep 
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track of holes in this sequence of cubical complexes, as they grow from complex (1) to complex (5). In the 
simple example of Figure persistent homology will allow us to deduce that a single one-dimensional hole 
can be observed starting with the third complex, and that it will disappear again at level (5). 

The example involving gray-scale images can easily be generalized. Formally, a filtration of cubical com¬ 
plexes is a nested sequence of cubical complexes 

/Cl C a :2 C ... C /Cm . ( 6 ) 

Notice that since we require every AC/ to be a cubical complex, if a cube of dimension k appears for the first 
time in AC/, all of its boundary elements have to appear at the latest in AC/, but they could appear earlier. 
While filtrations could be specified explicitly, one often convenient, and in fact equivalent, way of defining 

them is based on the concept of a filtering function. This is a function f : ACm —> {1. m} which assigns 

every cube C in the final cubical complex an integer between 1 and m in such a way that f(C) = / if and 
only if C G fCi and C 0 /C/_i, where we define ACo = 0. One can readily see that for this filtering function 
we have 

/C/ = {C G /Cm : f{C) < /} for all I < i < m . (7) 

Equivalently, but more useful for practical applications, assume that AC denotes a cubical complex and that 

we are given a function f : /C —>• {1. m} such that the following holds: If /\, 6 G AC are any two cubes 

such that B is contained in the boundary of A, then the function f satisfies f{A) > f{B). In this case, we 
call f a filtering function on AC, and Q defines a filtration of cubical complexes with /Cm = /C. We would 
like to point out that the above monotonicity property of f from a cube to its boundary is all that is needed 
to ensure that every AC/ is indeed a cubical complex. 

The above concept of a filtering function can easily be related to our image example from Figure In 
this case, let AC denote the cubical complex consisting of all sixteen pixels, together with their edges and 
vertices. Furthermore, suppose that the gray-scale intensities are given by the five values ai < 32 < ... < 35 . 
Then for every two-dimensional pixel /\ in /C we define f{A) = i, if its gray-scale intensity is given by 3 /. For 
an edge 6 G /C, we let f{B) be the smallest value of f on a pixel which contains B as an edge, and for a 
vertex C G AC, the value f{C) is the smallest value of f on an edge which contains C. Then Q leads to 

the filtration shown in Figure |^cj. This example also shows that choosing {1. m} as the index set for 

a filtration of cubical complexes can be done without loss of generality — through a simple transformation, 
any finite and discrete subset of M can be mapped onto this set. Alternatively, we could keep the values of 

the filtering function in the set {ai. 35 }, and index the cubical complexes by these distinct real numbers. 

This extension will be useful later on. 

Having introduced the concept of a filtration, we now turn our attention to persistent homology. Suppose 
we are given a filtration of cubical complexes as in Then for each dimension k and each complex in the 
sequence we obtain a homology group Hk{)Ci). Note that these homology groups are in general not ordered 
by inclusion, despite the fact that the underlying cubical complexes are. However, if c denotes an arbitrary 
/c-cycle in /C/, then one can readily see that c is also a /c-cycle in /C/+i. This implies that we can define a 
map from Hi^{JCi) to /-//<(/C/+i) by mapping the equivalence class of c in the first homology group to the 
equivalence class of c in the second one. In order words, the inclusions in Q lead to a sequence of maps 


/-//((/C2) ^ Hk{}Cm) , 

and these maps are in fact group homomorphisms. 

It was shown in the last section that if 7 G /7/((/C/) is a nontrivial homology class, and if c is a k-chain 
representing this class, then c encloses a /(-dimensional hole in /C/. Now consider the class 7 induced by c 
in /-//(^(/C/_|_i). If this equivalence class is nontrivial, then c still encloses a hole in /C/_|_i. On the other hand, 
if the equivalence class 7 is trivial, then the hole has been filled in by cubes that were added to /C/ during 
the formation of /C/_|_i. In other words, the homomorphisms induced by inclusion allow us to detect the 
disappearance of holes. Note, however, that holes can also disappear in a second way. Rather than being 
mapped to a trivial homology class, it is also possible that two different homology classes 71,72 G /7/((/C/) 
are mapped to the same nontrivial class in /7/((/C/+i), and this case will be discussed in more detail later. 
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Figure 4: Illustration of persistent homology. The six images show a filtration of cubical complexes in the 
plane, which leads to the following homology classes and lifespans. ( 1 ) In the complex /Ci, two connected 
components are born, i.e., two nontrivial zero-dimensional homology classes. (2) In /C 2 , another two con¬ 
nected components are born. (3) In the complex /C 3 , the two connected components born in IC 2 die, since 
each of them merges with a component which was born in /Ci — these latter two components still live 
on. (4) In /C 4 , however, the remaining two components merge, so one dies and the other survives. Since 
both were born at the same time, one can randomly choose which one survives. (5) In the complex /C 5 a 
1-dimensional cycle is created, which gives rise to a nontrivial homology class. ( 6 ) This 1-dimensional class 
dies in the final complex JCq. 


What about the appearance of holes? Suppose again that 7 6 Hk{}Ci) is a nontrivial homology class 
representing a /c-dimensional hole in /C,. One can easily see that if this hole was already present in /C/_i, 
then 7 has to be in the image of the inclusion-induced homomorphism from Hk{ICi-i) to Hk{JCi). In other 
words, the hole 7 appears for the first time in /C/, if and only if it is not in the range of this homomorphism. 

Being able to track the appearance and disappearance of holes in the filtration is the main idea behind 
persistent homology. While it is usually defined via the so-called persistent homology groups, which are just 
the images of the above-mentioned maps between homology groups induced by inclusion, we proceed right 
away to an equivalent formulation which is based on the lifespan of homology classes. We say that a k- 
dimensional homology class a is born in the complex ICj, if we have both a G Hk{ICj) and a ^ i{Hk{ICj-i)), 
where l : Hk{ICj-i) —>• Hk{JCj) is induced by inclusion. Moreover, we say that a k-dimensional homology 
class a dies in the complex /C^, if Hk{IC(,) is the first homology group where a either became trivial, or where 
it became identical to another class which was born earlier in the filtration. Note that the first of these 
cases, corresponds to the disappearance of a hole discussed above, while the second case addresses the issue 
of two nontrivial homology classes merging — in this case, the class that was created later in the filtration 
dies, while the one which was created earlier survives. In other words, merging effects can be treated in a 
unique way, through a process called the elder rule. 

Using the definitions of the last paragraph, persistent homology provides us with a list of nontrivial and 
independent homology classes, together with their birth times and death times. If 7 is one of these classes 
and if its birth and death times are given by j and I, respectively, then the interval [/', 1) is called the persistence 
interval associated with 7 , and the integer i — J is called its lifespan. If there are classes which are born, 
but never die, then we set i = 00 . In typical applications, persistent homology classes with relatively long 
lifespan are considered to be more important than the ones that have shorter lifespan. 

The lifespans of homology classes are illustrated in Figurej^ As described in the caption, this figure shows 
a filtration of six cubical complexes in the plane. Persistent homology yields five different homology classes. 
Four of these are in dimension zero, and their persistence intervals are given by [l,oo), [1,4), [2,3), and 
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Figure 5: Persistence diagrams for the filtration shown in Figure The left image shows the persistence 
diagram in dimension zero, which contains the points (l,oo) and (1,4) once, and the point (2,3) twice. 
These multiplicities are indicated by the numbers adjacent to the points. The right image contains the 
diagram for dimension one, which only contains the single point (5,6). 


again [2,3), respectively. The one nontrivial class in dimension one has persistence interval [5,6). 

The persistence intervals provide information about when homology classes are born, evolve, and die in 
the filtration. But they also can be used to recover the homologies of the involved cubical complexes. It can 
be shown that the Betti number of the complex JC; in the filtration is given by the number of persistence 
intervals in dimension k which contain /'. For example, in the situation of Figure|5 the complex /C 3 has Betti 
number/ 3 o(/C 3 ) = 2 , since only the intervals [l,oo) and [1,4) contain the index 3. We would like to pint 
out again, however, that persistent homology provides considerably more information than the sequence of 
Betti numbers of the complexes /C/, as it allows us to track for how long specific homology classes survive. 

While the persistence intervals encode all of the information provided by persistent homology, dealing with 
the intervals directly is fairly cumbersome. In practice, these intervals are visualized in a persistence diagram 
as follows. Fix a dimension k. Then the persistence diagram in dimension k \s a finite collection of points 
(possibly of multiplicity larger than one) in the plane such that (a, b) € is a point in the k^^ diagram if and 
only if persistent homology provides a nontrivial /c-dimensional homology class with interval [a,b). In other 
words, a persistence diagram is a multiset of points in which corresponds to the persistence intervals. 
Notice that all of these points lie above the diagonal in the first quadrant of For the filtration in Figure]^ 
the persistence diagrams in dimensions zero and one are shown in Figure [^ 

For our special situation of a filtration indexed by integers as in ([^, the points in a persistence diagram 
always have integer coordinates. However, if one considers a cubical complex /C together with a filtering 
function f : /C —)■ {ai. am}, one can consider the associated filtration 

/Cai c ICa 2 c ... c /Can, defined via /Ca, = {C € /C : f(C) < a,} . ( 8 ) 

In this more general setting, persistent homology produces persistence diagrams in a completely analogous 
way, but now with points of the form {aj, aft) G R^ which lie above the main diagonal. As mentioned before, 
the values aj could be the gray-scale intensities in an underlying image. What happens if these values are 
slightly perturbed? One would expect that new persistence diagram is only a slight perturbation of the old 
one, provided the perturbations are small enough. 

One of the main features of persistent homology is that this perturbation idea can be quantified in a most 
satisfying way. In order to do this, one needs to introduce a metric on the space of persistence diagrams, 
and there are a number of ways for doing this (TH]. For the purposes of this paper, we focus on the so-called 
bottleneck distance, which is closely related to the earth mover's distance discussed in [23l. To define this 
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Figure 6: For the persistence diagrams Di and D 2 shown in the first two images, the third image contains 
a matching : Di ^ D 2 in which points have to be moved only over fairly short distances. 


distance, consider two persistence diagrams Di and D 2 , both in the same dimension. Note, however, that 
these persistence diagrams could be created from two completely different filtrations. Let us now assume 
that in addition to the finitely many points corresponding to persistence intervals, we add all points of the 
form (x,x) with x G M to both diagrams, each of these points with infinite multiplicity. While this might 
seem strange at first glance, it is necessary for comparing persistence diagrams with different number of 
points. We denote the resulting new diagrams by Di and D 2 , respectively. After these preparations, we can 
finally define the bottleneck distance of the diagrams Di and D 2 , which will be denoted by Woo{Di, D 2 ) and 
is defined as 

l/l/oo(Di, D2) = inf „ ( sup ||x-ib(x)|| 

bijections b\bi^b2 \xeDi 

In other words, the bottleneck distance tries to find a matching b : Di —)■ D 2 which minimizes the largest 
distance between a point x G Di and its image b{x) G D 2 . Note that obtaining matchings of this type 
is made possible only through extending the persistence diagrams by the points on the diagonals, since Di 
and D 2 could have very different numbers of points. Since points on the diagonal would correspond to 
homology classes which die immediately when born, this should not alter the information content of the 
persistence diagrams. Intuitively, the bottleneck distance between two diagrams Di and D 2 is of size e, if it 
is possible to move every point of Di to a point in D 2 or the diagonal by less than distance e, as well as vice 
versa. This is illustrated in Figure As we mentioned before, typically, persistent homology classes with 
relatively long lifespan are considered to be more important than the ones that have shorter lifespan — and 
the latter ones give rise to points in the persistence diagram which are closer to the diagonal. This implies 
that these are the points that can easily be moved onto the diagonal in a matching. 

We now return to filtrations which are generated through a filtering function f : /C —>• {ai. am} as 

defined in How do perturbations of f affect the resulting persistence diagrams? To describe this, let 
us consider two filtering functions fi and f 2 on the same cubical complex /C, but we do not assume that 
the finite ranges of fi and f 2 are the same. If we further denote the persistence diagrams in dimension k 
associated with fi and by Dk{JC, fi) and Dk{JC, f 2 ), respectively, then one can prove a stability theorem in 
the form of the estimate 

W^{Dk{K:.fi).Dk{IC.f2)) < \\f 1 -f 2 \L . (9) 

where ||fi — /^||oo denotes the usual maximum norm. See [18] for more details. In other words, the bottleneck 
distance of the two persistence diagrams is bounded by the maximal function value difference of the two 
filtering functions. 

The stability theorem encapsulated by @ is one of the central properties of persistence diagrams, which 
is invaluable in applications. It implies, for example, that if observational data is perturbed by noise which is 
bounded by a constant e, then the persistence diagram of the noisy data will be at most e away from the 
persistence diagram of a true underlying data. In other words, persistence is a robust topological metric. 
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In the context of this paper, this robustness property allows us to compute persistence diagrams of finite¬ 
dimensional spectral approximations to the true solution u of the Cahn-Hilliard-Cook equation Q, which are 
basically indistinguishable from the true underlying persistence diagrams due to spectral accuracy. 

In closing, and for the sake of completeness, we would like to point out that the bottleneck distance is 
only one possibility for measuring the distance between persistence diagrams. In fact, in some situations its 
insistence on only registering the largest move distance ||x — b(x)|| is too stringent. In such cases one can 
make use of the Wasserstein distance HE], which is defined as 


Wp{Di,D2) = ^ ik-K^)ri 

bijections b-.Di^D2 \ ~ / 

\xgDi / 

for p > 1 . This distance has the advantage of taking all line segments in a matching into account, but it 
comes at the price of no longer satisfying the stability result However, more elaborate stability estimates 
which are valid for the Wasserstein distance are available in a number of situations, and we refer the reader 
to m for more information. 

2.3 Persistence Landscapes 

By moving from the concept of homology groups to the concept of persistent homology, we have made 
the transition from thresholding a gray-scale image at one specific level to considering the fine structure of 
the image as we pass through all possible shades of gray. As we mentioned in the introduction, the goal of 
this paper is to show that the topology of material microstructures such as the ones in Figure encodes 
interesting facts about the underlying dynamics. After suitable discretization these microstructures can be 
represented as gray-scale images, with the value of the solution u of @ serving as the filtration parameter. 
Recall, however, that the Cahn-Flilliard-Cook model is a stochastic partial differential equation, and therefore 
with probability one no two realizations of the dynamics will lead to identical patterns. In order to uncover 
this “typical" behavior, one would therefore like to consider statistical properties, in particular, averages of 
the topological information. 

While averages of Betti numbers are easily accessible, and have in fact been used for example in [T7], it is 
not immediately clear how the average of persistence diagrams can be computed. One of the first attempts 
at defining the mean of two persistence diagrams is the notion of Frechet mean which was explored in EH 
Without going into too much detail, assume we are given two persistence diagrams Di and D 2 , and assume 
that we have found a matching b : Di —)■ D 2 which realizes the p-Wasserstein distance Wp{Di, D 2 ), for 
some p > 1. Then the Frechet mean of the diagrams Di and D 2 is defined as the collection of midpoints 
of the line segments between x and b(x), for all x G Di. Despite the fact that the Wasserstein distance 
between Di and D 2 takes the length of all matching distances ||x — b(x)|| into account, one can easily see 
that the above notion of Frechet mean is not well-defined. Consider for example a diagram Di which consists 
of the points (J, k) and 0 + 1, k + 1), while the diagram D 2 has the points {j, k + 1) and (J + 1, k), and where 
we assume that j+1 < k. Then there are two minimal matchings which involve these points: One matching 
maps along two horizontal line segments of length one, while the other one uses two vertical line segments 
of length one. This implies that in the former case the Frechet mean consists of the points {J + 1/2, k) 
and 0 + 1/2, k + 1), whereas in the latter case it is given by the points (J, k + 1/2) and 0 + 1. k + 1/2). 

It is not surprising that this obvious non-uniqueness of the Frechet mean makes a statistical analysis 
of persistence diagrams which involves this concept nontrivial. From a mathematical point of view, these 
problems can be resolved, see for example m Unfortunately, however, these ideas have not yet been 
implemented algorithmically, as they pose serious technical and computational difficulties. For the purposes of 
the present paper we therefore pursue a different approach, which uses the concept of persistence landscapes. 
These objects were introduced in [5], and they prove to be an algorithmically feasible tool for the statistical 
analysis of persistence diagrams. 

In order to motivate the definition of persistence landscapes, we return one final time to the Frechet 
mean. Its non-uniqueness stems from the fact that it operates on a discrete set of points. In some sense. 
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this situation strongly resembles the problem of computing the mean of two integers within the ring of 
integers. Also in this case, the mean is often not unique. In the integer case, however, the mean problem 
can be solved easily by introducing rational numbers. We will see below that the viewpoint behind persistence 
landscapes is very similar — it is based on the idea of embedding persistence diagrams into the larger space of 
piecewise linear functions, where the mean is uniquely defined. For the remainder of this section, we present 
an intuitive introduction to persistence landscapes, rather than a formal definition. For a more detailed 
discussion we refer the reader to [5]. Important for our applications is the fact that persistence landscapes 
can be computed and manipulated efficiently, and algorithms for this can be found in | 6 ]. 

To begin with, suppose we are given a point {b, d) G in a persistence diagram, i.e., the point corresponds 
to a homology class with birth time b and death time d. With this point, we associate a piecewise linear 
function f{b,d) : K ^ [0, oo), which is defined as 




0 

if 

x0(b, d) , 

x-b 

if 

X G {b. 

d-x 

if 

X G (^, d) 


( 10 ) 


In other words, the function f^b^d) is a hat function, whose maximal value is {d — b)/2, i.e., half the lifespan 
of the homology class represented by the point {b, d). Moreover, this global maximum of f(^b,d) occurs at the 
center of the associated persistence interval, and on either side the function decays linearly towards zero with 
slopes ±1. In fact, there is another, more geometric way of describing f(^b.d)- For this, map the point {b, d) 
in the persistence diagram via the linear map 


( b\ ( 1/2 1/2 

W J V -1/2 1/2 



( 11 ) 


and then draw lines with slopes ±1 from the image of the point towards the horizontal axis. Clearly, the 
function f(^b,d) encodes the original data provided by the point {b,d), but in a more convenient way: The 
support of the function equals the closure of the persistence interval [b, d), while the height of the function 
is half the lifespan. Note also that the map ( fTI| ) transforms the main diagonal into the horizontal axis. This 
procedure is illustrated in Figure [^a/, (b), and (c). 

After these preparations, we can now introduce the concept of persistence landscape. As was mentioned 
earlier, in many applications a homology class is considered important, if its lifespan d — b \s large. Since 
the heights of the functions f(b,d) shown in Figure [^c/ are proportional to the lifespan, the functions close 
to the top are therefore considered more significant than the ones closer to the horizontal axis. Motivated 

by this observation, the persistence landscape of the birth-death pairs (bi.d;), where / = 1 . m, which 

constitute the given persistence diagram is the sequence of functions A/< : M —)■ [0, oo) for /c € N, where Xk{x) 

denotes the largest value of the numbers f^bi.dpi^)’ for / = 1 . m, and we define Xk{x) = 0 \f k > m. 

Equivalently, this sequence of functions can be combined into a single function : N x M — )■ [0, oo) of two 
variables, if we define L{k, t) = Xk{t). See again Figure|^for an illustration of these concepts. 

One can easily see that the persistence landscape L encodes exactly the same information as the underlying 
persistence diagram. Yet, in some sense, this information is presented in a way which is mathematically more 
accessible. First of all, it creates a filtered representation of the persistence diagram which is premised on the 
assumption that homology classes with large lifespans are more important than ones with shorter lifespans. 
This is reflected in part by the monotonicity property 


L{k.x) > L{k + l.x) 


for all 


/( e N and x G 


( 12 ) 


Secondly, and most importantly, it represents the persistence data as a piecewise linear function L, and 
averaging functions of this type results again in a piecewise linear function. Moreover, the monotonicity 
property ( fl^ is transferred from the involved landscapes to the average. 

Thus, the introduction of persistence landscapes allows us to easily talk about the average of persistence 
information. Assume we are given persistence landscapes Li . L^. Then the averaged persistence land¬ 

scape is defined by X)/=i l-i/£, and it is again a piecewise linear function which satisfies (12). In addition. 
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Figure 7: The idea behind persistence landscapes: (a) contains the initial persistence diagram, while (b) shows 
the transformed diagram after mapping all points and the diagonal with the transformation (c) depicts 
the piecewise linear functions f(b.d) defined in (10) for every one of the transformed points. Finally, the graph 
in (d) illustrates the first persistence landscape function Ai in red, which is the upper envelope of the graphs 
from the previous picture. The second and third persistence landscape functions A 2 and A 3 are shown in 
green and blue in (e) and (f), respectively, and they follow the second- and third-largest function values of 
the functions f(^b,d) image (c). All subsequent persistence landscape functions A/< are identically zero. 


it is clearly uniquely defined. But considering the persistence information in the form of a real-valued func¬ 
tion defined on N x E has another advantage. For such functions, we can use the standard /.'^-metric to 
quantify how “different” two landscape functions are. In fact, it was shown in [5] that persistence landscapes 
are stable with respect to the /.'^-metric, i.e., small changes in the filtering function which generates the 
underlying persistence diagram leads to small changes in the resulting persistence landscape with respect to 
the /.'^-norm. 

Also from a computational point of view persistence landscapes are useful. The construction of persistence 
landscapes has been implemented by the first author of this paper in a library called Persistence Landscape 
Toolbox, see [S] for more details. The library is publicly available and can be downloaded from the first 
author’s webpage. It implements various manipulations of persistence landscapes, including the statistical 
classifier which will be heavily used in the analysis of the present paper, and which will be described in more 
detail in the next section. 


3 Topological Analysis of the Cahn-Hilliard-Cook model 

We now turn to the analysis of the Cahn-Flilliard-Cook model using persistence landscapes. We begin in 
Section |3.1| by introducing the basic methodology, including a description of the persistence landscapes 
based statistical classifier and of the numerical simulation techniques. After that. Section is concerned 


illustrates that even the decomposition stage can be recovered accurately from the topological information. 


with showing that the topology evolution of (Q encodes the central mass parameter n, while Section 33 
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3.1 Basic Methodology 

In order to generate the microstructures for the persistence landscape analysis, we simulated the Cahn- 
Hilliard-Cook model @ at the parameter values 

e = 0.005 , a = 0.001 , with space-time white noise £ , 

and on the square domain Q = (0, Vf- c All simulations were performed in C++ using a linearly implicit 
spectral method. More precisely, the solution u{t,x) for f > 0 and x = (xi,X 2 ) e Q is approximated using 
the Fourier expansion 


K-l 

u{t,x) « ^ Uki.k 2 {t) ■ Ql.te • COS/CiTTXi • COS/C27rX2 

/Ci,/C2=0 


(13) 


where the constants chosen in such a way that the functions cos /cittxi cos/(27rx2 form a 

complete orthonormal set in for ki, k 2 € Nq. Also the white noise process £ is approximated in this 

way, i.e., we use cut-off noise which acts independently on all K^ — 1 Fourier modes in the representation (13). 


Recall that the noise is assumed to be mass-conserving, so it does not act on the constant mode. Throughout 
this paper, the simulations use K = 256, and to avoid aliasing effects, the necessary Fourier coefficients 
of the nonlinearity F'{u) are computed using the two-dimensional discrete cosine transform on a grid of 
size 512^. 

For a variety of different mass values n, we then performed the following simulations. Let uq g L^{^) 
denote a randomly chosen initial condition with average n which satisfies 


= 10 


Then the above numerical method is used to compute the solution u of the Cahn-Hilliard-Cook model 
with iv(0, •) = Uq for times 


0 < f < Te 


80e2 

FW 


@ 


where we recall that F"{ii) = 3y?- — 1. While the choice of simulation endtime might seem strange at first 
sight, it takes into account the different instability strengths of the homogeneous state ii as a function of 
the mass. In fact, by choosing the endtime Tg as above every simulation leads to basically the same phase 
separation horizon. The time interval is discretized using 100000 steps, and solution snapshots are recorded 
at a variety of times, as will be described in more detail below. 

The solution snapshots u{U, ■) obtained through this procedure are then analyzed using persistence land¬ 
scapes. As a first step, the spatial domain is discretized into 512^ equal-size subsquares, and on each square 
the function u is approximated using its function value at the center of the subsquare. In this way, we obtain 
a pixelated image version of the microstructure, as shown in Figure [T] where the function values of u take 
the role of the gray-scales in the previous section. Next, the solution range [—1,1] is discretized into 256 
equal-sized intervals, and the persistence diagrams are computed by considering sublevel sets of u{U, •) as 
the threshold parameter increases in 256 steps from —1 to 1. Notice that since the solution u could take 
values slightly larger than one, we actually consider the function minju, 1} instead of u. In the language 
of Section |2.2| we consider the cubical complex fC which consists of 512^ equal squares to form the unit 
square Q (where one unit in the cubical complex representation corresponds to 1/512 units in the Q-scale), 
and the filtering function min{u, 1} : Q —>• (—oo, 1] evaluated at the centers of the subsquares. Together 
with Q this leads to a filtration /Ci C ... C /C 256 , where the thresholds a,- are chosen uniformly between —1 
and -Fl. The persistence diagrams are computed using the software PHAT [T], together with an implemen¬ 
tation of cubical complexes. Those computations can also be performed with Perseus m The associated 
persistence landscapes are then determined using I^. Some of the resulting persistence landscapes can be 
found in Figures and for dimensions zero and one, respectively. As we pointed out in the last section, 
none of the patterns discussed in our setting will have non-trivial two-dimensional homology, as the underlying 
two-dimensional cubical complexes are flat. 
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Figure 8: Sample persistence landscapes in dimension zero for microstructures created by the Cahn-Hilliard- 
Cook model From top to bottom, the rows correspond to mass ^ = 0, 0.01, and 0.02, respectively. 
From left to right, the respective columns use solution snapshots at times t = 0.2Te, O.QTq, and Tq. 



Figure 9: Sample persistence landscapes in dimension one for microstructures created by the Cahn-Hilliard- 
Cook model From top to bottom, the rows correspond to mass jU = 0, 0.01, and 0.02, respectively. 
From left to right, the respective columns use solution snapshots at times t = 0.27’e, O.OTe, and Tg. 
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Due to the stochastic nature of the Cahn-Hilliard-Cook model, the above simulations are repeated many 
times in a Monte-Carlo type fashion. In this way, we created three data sets that will then be analyzed 
further: 


(Dl) The first data set considers the 51 mass values jj. = 0.01 r? for n = 0,..., 50, and in each case five 
different solution paths of Q, i.e., five different initial conditions uq. Along the solution, snapshots 
are saved at times t = 0.002iTe for £ = 1 ,..., 500. 

(D 2 ) The second data set considers the 6 mass values jU, = 0.1 r? for r? = 0,..., 5, and in each case fifty 
different solution paths of Along all solution paths, snapshots are saved at times t = 0.002£Te 
for £= 1.500. 

(D3) The third and last data set considers the 51 mass values ijl = 0.01 r? for n = 0,..., 50, and in each 
case 100 different solution paths of Along each solution, snapshots are saved at times t = 0.04£Te 
for £= 1.25. 


Thus, in each of the three data sets we consider N different mass values jU, where N = 51,6,51, respec¬ 
tively, and for each simulation we obtain solution snapshots at M different times, where M = 500,500,25, 
respectively. Every one of these solution snapshots leads to two associated persistence landscapes, one each 
in dimensions zero and one. 

We would like to briefly comment on the fact that all of our data sets consider non-negative mass 
values M > 0. Based on our thresholding process described above, we study sublevel sets of the function u 
as the threshold increases from —1 to +1. In Figure[^ these sublevel sets for the threshold 0 = 0 are shown 
in red. Thus, for mass values 0 jj, < 0.5 the sublevel sets always have nontrivial 1-dimensional homology, 
which in fact measures the bulk or interior behavior of the material, as described in [IlllTT]. I n contrast, 
one can easily see that for mass values -0.5 < < 0 the sublevel sets usually have trivial first homology, 

i.e., the first homology dimension is useless for topological classification. 

Based on the above topological data, our goal in the next two sections is a statistical study to solve inverse 
problems for the total mass ii and the decomposition stage. While this could be done in a number of ways, we 
concentrate on one statistical classification method which has been implemented in and which we recall 

now. Suppose we are given N different classes Ci.C/v of persistence landscapes. Assume further that 

each class Cn consists of S different landscapes Lnj, for j = 1 .S. More precisely, in every class €„ we 

have S different persistence landscapes L^j, j = 1.S, in dimension k, where k 6 Nq. Our classification 

scheme will be based on the previously mentioned /.'^-distances between persistence landscapes, and therefore 
we choose and fix a real number p G [ 1 , oo]. Suppose we are given a new persistence landscape L and would 

like to decide which of the classes Q.C/v it belongs to, based only on the persistence information in 

dimension k G Nq. Then we use the following classification scheme: 

(Ck) Suppose N classes of persistence landscapes are given as above. Then for each n = 1 .A/ we define 

the average classifier of the class in dimension k via 


c' = 1 • V . 

7=1 


(14) 


Now let L be any other persistence landscape which we would like to classify. Then we say that L has 
been classified to belong to class Cn using dimension k, if there exists a unique index r? G {1,..., A/} 
such that 


^ n 


< 


LP 




LP 


for all 


J = l- 


N 


(15) 


If no such unique index n exists, then we say that the classifier fails to classify L. 


Intuitively, the classification scheme (Ck) is easy to understand. Based on the different persistence landscapes 
available in a given class, we compute their average to determine the “typical" persistence behavior in this 


18 










class. Given a new persistence landscape L, one then uses the average classifier which lies closest to L to 
determine the class which most likely contains L. 

The above classification scheme (Ck) in fact describes a whole family of classification schemes, indexed 
by the homology dimension k. In our situation, only the cases k = 0 and k = 1 lead to nontrivial persistence 
landscapes, and therefore we will use only schemes (CO) and (Cl) in the following. Note, however, that in 
some situations one might like to use all available topological information, i.e., all dimensions k, at the same 
time to determine the most likely class which contains a given persistence landscape L. This can be achieved 
by the following simple scheme. 




(CA) Suppose N classes of persistence landscapes are given as above, and that the average classifier of 


the class in dimension k has been defined as in ( |14D , for r? = 1 ,..., A/ and k G Nq. As before, let L 
be any other persistence landscape which we would like to classify. 


For each dimension k G No, one can then order the N possible classification classes according to 
monotone increasing values of the distances ||/.''— C„||/.p, and we let Sj denote the first J class indices 

in this sequence, for j = 1. N. Note that is a singleton set which contains a class whose 

averaged classifier C„ is closest to , while we have 5^ = {1. N} for all k G Nq. Furthermore, 

let ^ = 1. N be the smallest number such that 0. 

Then we say that L has been classified to belong to class €„ using all dimensions, if any ordering of 
the classes as above leads to 

OO 

n = {^} ■ (16) 

k=0 

If this intersection contains more than one element, or if there are ties in the above orderings due to 
equidistant average classifiers, then we say that the classifier fails to classify L. 


This classification scheme can be useful if no clear choice of dimension k in (Ck) is evident a-priori, and we 
will use it also in our further studies. 

For our application to the Cahn-Hilliard-Cook equation in Section we are using a slight extension 
of the above classification schemes. For each simulation, our input data consists of a sequence of M 

different persistence landscapes which correspond to microstructure patterns at times ti . tM- Such a 

sequence of persistence landscapes will be called a topological process. If we define the average of topological 

processes Pi. Ps as the topological process whose persistence landscape in dimension k is the average 

of all the persistence landscapes in dimension k of the processes Pi.Ps, then we can easily talk 

about an averaged classifier as in (Ck). Furthermore, if we define the /.'^-distance between topological 
processes P and Q as the sum of the /.^-distances between the persistence landscapes in the F'^ positions of 

the processes P and Q, for / = 1. M, then one can easily reformulate for the case of topological 

processes. The total classification scheme (CA) can be extended analogously. The resulting schemes will be 
used for the classification of the Cahn-Hilliard-Cook simulations described in (Dl), (D2), and (D3). 


3.2 Identifying the Mass via Topology 

As described in the last section, our simulations of the Cahn-Hilliard-Cook model © provide snapshots of a 

solution realization along M equidistant timesteps ti. tM- For each time t^, computational topology is 

then used to create the associated persistence landscapes in dimensions 0 and 1. In other words, for each 
simulation of © we obtain a topological process which captures the persistence evolution of the created 
microstructure over the time interval (0,Te]. In this sense, our topological data is a refinement of the 
topology evolution curves which were considered in m 

What does this refined topology evolution capture? In the present section, we demonstrate that it is in 
fact possible to determine the total mass p, of the underlying simulation extremely accurately. This implies 
the somewhat surprising observation that the total mass p has quite a significant influence on the pattern 
formation process — to the extent that persistence information alone, which of course only depends on 
/(-dimensional connectivity measurements, but not on size information, suffices to accurately infer p. 
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L^-norm 

Batch No. 

Flits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Wrong 

1 

123 

29 

0 

0 

1 

0 

2 

128 

23 

1 

0 

1 

0 

/.^-norm 

Batch No. 

Flits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Wrong 

1 

112 

39 

1 

0 

1 

0 

2 

110 

41 

1 

0 

1 

0 

L^-norm 

Batch No. 

Flits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Wrong 

1 

78 

58 

14 

1 

1 

1 

2 

82 

58 

12 

0 

0 

1 


Table 1 : Classification results for the data set (Dl) using the classification scheme (CA). The two training 

_ ^ 

set batches consist of R 7 = 2 simulations each. Once the averaged classifier processes for /? = 1 . N, 

where N = 51, have been computed, the remaining {R — Rt)N = 153 solution snapshot processes are 
classified using the L^-, L^-, and /.°°-norms. The table contains the number of classification hits as well as 
detailed information on the index difference observed for the misclassifications. The column labeled “Wrong" 
contains classifications which missed by more than 4 indices. 


We demonstrate this mass-topology dependence using the three data sets (Dl), (D 2 ), and (D3) described 
The simulations in each of these data sets give rise to topological processes of length M, and 


in Section 3.1 


they can be divided into N different mass classes. While in (Dl) and (D3) we considered the A/ = 51 mass 

values 11 = 0.01 /? for /? = 0.50, the data set (D2) considers the A/ = 6 mass values 11 = 0.1 n, where now 

we have n = 0 .5. Furthermore, the length of the topological processes is given by M = 500 for (Dl) 

and (D2), and M = 25 for data set (D3). For each mass value n, a data set contains R stochastically 
independent repetitions of the simulation. For (Dl), (D2), and (D3) we have R = 5, R = 50, and R = 100, 
respectively. 

Based on the composition of the data sets, we apply the classification schemes (CA), as well as (CO) 

and (Cl), from the previous section for mass value classes Ci.C/v, which correspond to the chosen N 

mass values iJ, in increasing order. Then we proceed as follows. 

• Based on Rt < R training runs each from the N considered mass values, we determine the averaged 
classifier of the class in dimension k as defined in (14) with S = Rj, but extended for the case 
of topological processes. 


• We then try to classify each of the remaining (/? — Rt)R simulations in the data set, by either using 
the classification scheme (CA), or one of the schemes (CO) or (Cl). 

In each case, we note the index n obtained from the classification scheme, and compare it to the actual mass 
index ritrue of the underlying simulation. If /? = ritrue. then we have obtained a "Flit”, otherwise we say that 
the classification “missed by \n-ntrue\"- Note that the integer value \n — ntrue\ is directly proportional to the 
mass difference incurred by the classification. For the data sets (Dl) and (D3) the associated proportionality 
factor is 0.01, while for (D2) it is 0.1. 

We now turn to the first data set (Dl), which consists of R = 5 simulations each for N = 51 different 
mass values. In this case, we took the first Rj = 2 data sets for each mass value as training sets, and 
their average constitutes the averaged classifier in each case. Classifying the remaining (/? — Rt)R = 153 
topological process using scheme (CA) as described above then leads to the results shown in Tablein the 
rows denoted “Batch 1". If instead we use the third and fourth simulations as two-element training set, one 
obtains the lines labeled as "Batch 2". For the classifications, we used the /.'’-norm with p = 1, 2, 00 . 

The results in the table are remarkable. If one uses the /.^-norm for the classification, then in over 80% 
of the cases the classification scheme finds the correct mass value class. If one allows for an error of the 
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Dimension k = 0 

Norm 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Missed by 5 

Wrong 

U 

112 

40 

0 

0 

0 

0 

1 

U 

100 

48 

4 

0 

0 

0 

1 

1^00 

81 

53 

11 

3 

1 

0 

4 

Dimension k = 1 

Norm 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Missed by 5 

Wrong 

U 

129 

22 

1 

0 

1 

0 

0 

U 

120 

29 

3 

0 

1 

0 

0 

1^00 

68 

59 

18 

4 

3 

1 

0 


Table 2: Classification results for the data set (Dl) using the classification schemes (CO) and (Cl). The 
training set consists of two simulations, the remaining three simulations for each of the A/ = 51 mass values 
are then classified using the L^-, L^-, and /.°°-norms. The table contains the number of classification hits as 
well as detailed information on the index difference observed for the misclassifications. The column labeled 
“Wrong" contains classifications which missed by more than 5 indices. 


form jLi ± 0.01, then the classification succeeds in at least 98% of the cases. While the behavior for the 
/.^-norm is only slightly worse, the /.°°-norm does not seem to be quite as accurate an estimator. The 
latter norm gives correct classifications only in at least 50% of the cases, errors of the form ii ± 0.01 in 
at least 88%, as well as errors n ± 0.02 in at least 98% of the cases. In addition, the two classifications 
appearing in the “Wrong" column have index differences of 21, i.e., they are significant misclassifications. 
Despite the worse performance of the /.°°-norm estimator, these results clearly indicate how powerful the 
above topological classification scheme is. We would like to point out that one cannot expect the classifier 
to have 100% accuracy. On the one hand, these classifications used training sets which consist of two 
topological processes, which clearly is an extremely small batch size. Moreover, we are classifying pattern 
evolutions created by a stochastic partial differential equation which originate at random initial conditions. 
Given these uncertainties, the presented topological classification method is highly effective in distinguishing 
the pattern evolutions from the Cahn-Hilliard-Cook equation 

For the above classifications we have used the classification scheme (CA) which simultaneously considers 
persistence information in dimensions k = 0 and k = 1. How do the dimension-dependent classification 
schemes (CO) and (Cl) fare? The results for "Batch 1” as above are collected in Table As before, both 
the L^- and the L^-norm classifications lead to precise mass value classifications, while the classification based 
on the /_^-norm is worse. Note, however, that the performance of the (Cl) scheme is consistently better 
than the one of the (CO) scheme. We believe that this can be explained as follows. In our simulations, we 
considered mass values ^ > 0. This implies that in the droplet forming regime 0 <C M < 0.5, most threshold 
values lead to a connected cubical complex with a number of holes, which correspond to the droplets of 
one material forming in the background matrix consisting of the second material. In this situation, the first 
Betti number only counts the interior droplets, i.e., the ones that do not touch the boundary. In other 
words, in our setup, the 1-dimensional Betti number of the considered sublevel set is equal to the number 
of interior components of the minority phase. In contrast, the 0-dimensional Betti number counts both 
the interior and the boundary components of the majority phase. Thus, it seems reasonable to attribute 
the better performance of the scheme (Cl) to its exclusive focus on the bulk behavior of the material, 
while (CO) measures the combination of both bulk and boundary phenomena. As was pointed out in |T4] . 
these phenomena exhibit different scaling behaviors, and can therefore affect the results for certain values 
of the parameter e > 0. 

Before we move on to the other two data sets, we would like to elaborate a bit more on the relative 

behavior of the /.^-classifiers for p = l,2,oo in the above setting. One would expect good classification 

_ 

results only if the averaged classifiers defined in (14), but extended to the case of topological processes, 

istance. More precisely, one would expect that the 


are suitably far apart with respect to the considered d 
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Figure 10: Distances between the averaged classifiers of the mass value class for the data set (Dl), 

where n = 1 .51. The horizontal and vertical axes correspond to the mass value index n, and the sum of 

the distances for dimensions k = 0 and k = 1, as defined in (]T), are color-coded as shown in the colorbars. 
From top left to bottom middle the images are for the L^-, L^-, and /.°°-norm, respectively. 


distances 


-C° 

^ n ^ m 


LP 


+ 


^ n ^ m 


LP 


for m, n = 1,... ,51 


(17) 


should be smoothly increasing as the distance \m — n\ becomes larger. The distances defined in ( [TtI ) for 
p = 1,2,00 are shown as heat maps in Figure [T^ from top left to bottom middle. Notice that while in 
all three cases this “monotonicity" can be observed, the level of smoothness of the heat maps decreases 
with increasing p. The first of these observations is responsible for the fact that all three choices of p lead 
to acceptable classification schemes, while the second observation explains the relative difference in their 
performances. 

We now turn our attention to the second data set (D2). In this case, we consider R = 50 simulations each 
for A/ = 6 different mass values, and the length of each topological process is still given by M = 500. We 
increase the size of the training sets to Rt = 10, and classify the remaining (/? — Rt)R = 240 simulations 
using the classification schemes (CA), (CO), and (Cl). Regardless of the choice of L'^-norm, where again we 
use p = 1, 2, oo, the classification schemes yields 100% accuracy. We have also performed cross validations 
by choosing different training sets of size ten, which again leads to perfect classifications. These results 
certainly owe to the fact that the discrete mass values p, considered in the data set are spaced further apart. 
However, while in the case (Dl) the /.°°-classification led to isolated misclassifications with large deviations, 
this was not observed in the data set (D2). We believe that these outliers are controlled due to the increased 
size of the training sets. 

The results so far have used fairly fine sampling in the temporal direction, which has led to two data sets 
with topological processes of length M = 500. It is natural to wonder whether the classification schemes 
work as well if this sampling size is decreased. To study this, we consider the data set (D3), which consists 
of R = 100 simulations each for N = 51 different mass values, but only saves M = 25 solution snapshots in 
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L^-norm 

Batch No. 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Wrong 

1 

2340 

210 

0 

0 

0 

0 

2 

2328 

222 

0 

0 

0 

0 

/.^-norm 

Batch No. 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Wrong 

1 

2188 

361 

1 

0 

0 

0 

2 

2163 

386 

1 

0 

0 

0 

L^-norm 

Batch No. 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Wrong 

1 

1655 

769 

113 

11 

2 

0 

2 

1699 

721 

117 

11 

2 

0 


Table 3: Classification results for the data set (D3) using the classification scheme (CA). The two training 

_ 

set batches consist of Rj = 50 simulations each. Once the averaged classifier processes C„ for /? = 1. N, 

where N = 51, have been computed, the remaining {R — Rt)N = 2550 solution snapshot processes are 
classified using the L^-, L^-, and /.°°-norms. The table contains the number of classification hits as well as 
detailed information on the index difference observed for the misclassifications. The column labeled “Wrong" 
contains classifications which missed by more than 4 indices. 


Dimension k = 0 

Norm 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Missed by 5 

Wrong 


2113 

433 

4 

0 

0 

0 

0 


1857 

669 

24 

0 

0 

0 

0 

1^00 

1581 

781 

142 

32 

7 

3 

4 

Dimension k = 1 

Norm 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Missed by 5 

Wrong 


2308 

242 

0 

0 

0 

0 

0 


2210 

315 

25 

0 

0 

0 

0 

1^00 

1547 

643 

235 

80 

28 

10 

7 


Table 4: Classification results for the data set (D3) using the classification schemes (CO) and (Cl). The 
training set consists of Rt = 50 simulations, the remaining 50 simulations for each of the N = 51 mass 
values are then classified using the L^-, L^-, and /.°°-norms. The table contains the number of classification 
hits as well as detailed information on the index difference observed for the misclassifications. The column 
labeled "Wrong" contains classifications which missed by more than 5 indices. 


the interval (0, Tg]. As a first experiment, we use two different training sets of size Rt = 50, and in each 
case classify the remaining (/? — Rt)N = 2550 simulations. The two different sets are referred to as “Batch 
No. 1" and "Batch No. 2”, respectively, the results for the classifier (CA) are collected in Tablewhile the 
results for (CO) and (Cl) can be found in Table 

For the classification scheme (CA), the simulation results are convincing. If one uses the /.^-norm for the 
classification, then in over 91% of the cases the classification scheme finds the correct mass value class, and 
all classifications are correct up to an error of the form iiztO.Ol. As before, the /.^-norm classifier performs 
slightly worse, with correct classifications in over 84% of the cases, and errors of ± 0.01 in over 99% of 
the cases. The maximally observed errors where of the form n ± 0.02. While the integral norms perform 
extremely well, the maximum norm /.°°-estimator again does not seem to be quite as accurate. It leads to 
correct classifications only in at least 64% of the cases, errors of the form fj, ± 0.01 in at least 94%, as well 
as errors jLi±0.02 in at least 99% of the cases. We would like to point out, however, that in all of the above 
situations, the performance exceeded the one observed in Table for the data set (Dl), despite the fact 
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L^-norm 

Batch No. 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Missed by 5 

Wrong 

1 

4174 

416 

0 

0 

0 

0 

0 

2 

4179 

411 

0 

0 

0 

0 

0 

3 

4163 

427 

0 

0 

0 

0 

0 

4 

4181 

409 

0 

0 

0 

0 

0 

5 

4143 

447 

0 

0 

0 

0 

0 

6 

4146 

444 

0 

0 

0 

0 

0 

7 

4149 

441 

0 

0 

0 

0 

0 

8 

4198 

392 

0 

0 

0 

0 

0 

9 

4173 

416 

1 

0 

0 

0 

0 

10 

4175 

415 

0 

0 

0 

0 

0 

L^-norm 

Batch No. 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Missed by 5 

Wrong 

1 

3859 

719 

12 

0 

0 

0 

0 

2 

3860 

724 

6 

0 

0 

0 

0 

3 

3809 

771 

10 

0 

0 

0 

0 

4 

3871 

713 

6 

0 

0 

0 

0 

5 

3747 

835 

8 

0 

0 

0 

0 

6 

3852 

734 

4 

0 

0 

0 

0 

7 

3846 

741 

3 

0 

0 

0 

0 

8 

3866 

721 

3 

0 

0 

0 

0 

9 

3843 

743 

4 

0 

0 

0 

0 

10 

3858 

729 

3 

0 

0 

0 

0 

/.°°-norm 

Batch No. 

Hits 

Missed by 1 

Missed by 2 

Missed by 3 

Missed by 4 

Missed by 5 

Wrong 

1 

3019 

1350 

205 

11 

4 

1 

0 

2 

3004 

1373 

185 

23 

5 

0 

0 

3 

2972 

1395 

207 

14 

2 

0 

0 

4 

2899 

1417 

255 

14 

5 

0 

0 

5 

2903 

1430 

236 

17 

4 

0 

0 

6 

2954 

1412 

195 

23 

5 

1 

0 

7 

2934 

1403 

234 

16 

3 

0 

0 

8 

2895 

1465 

213 

14 

3 

0 

0 

9 

2908 

1432 

231 

17 

2 

0 

0 

10 

2977 

1382 

203 

23 

4 

1 

0 


Table 5: Classification results for the data set (D3) using the classification scheme (CA). The ten training 

—k 

set batches consist of 10 simulations each. Once the averaged classifier processes for n = 1. N, 

where N = 51, have been computed, the remaining 90N = 4590 solution snapshot processes are classified 
using the L^-, L^-, and /.°°-norms. The table contains the number of classification hits as well as detailed 
information on the index difference observed for the misclassifications. The column labeled “Wrong" contains 
classifications which missed by more than 5 indices. 


that now our topological processes have only length M = 25. This increased accuracy is due to the larger 
training set size Rt = 50. 

How do the results change if we use the dimension-dependent classification schemes (CO) or (Cl)? The 
corresponding results are presented in Tablej^ and they largely confirm the observations made in the context 
of Table|^ Using only one dimension for the classification yields still acceptable, but definitely worse results. 
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While now the /.^-classification is clearly the most accurate, the L^- and /.°°-norm classifications come in 
second and third, respectively. Finally, the classification using dimension k = 1 performs better than the one 
in dimension zero, at least for the case of /.^-classifications. These results show, however, that in general it 
is preferable to use all available topological information for the classification. 

As a final experiment, we again consider data set (D3), but this time the R = 100 simulations for each 
mass value ii are partitioned into ten different training sets of size Rt = 10 each. These training sets are 

referred to as “Batch No. t for £ = 1.10, and for each training set the remaining {R — Rt)N = 4590 

topological processes are then classified using scheme (CA). The results of these computations are shown 
in Tableand they confirm our previous observations. For classifications using the L^-, L^-, and /.°°-norm 
one obtains exact hits in at least 90%, 81%, and 63% of the cases, maximal errors of the form ^ ± 0.01 
are observed in at least 99%, 99%, and 94% of the cases, and finally maximal errors ju ± 0.02 occur in at 
least 100%, 100%, and 99% of the cases, respectively. For the L^-, L^-, and /.°°-norm classifications, the 
maximal observed classification errors where ju,±0.02, ju,±0.02, and ju,±0.05, respectively. While all of these 
numbers are slightly decreased from the ones shown in Table this is to be expected due to the smaller 
training set size Rj = 10. Nevertheless, especially the /.^-norm classification performs amazingly accurate. 


3.3 Recovering the Snapshot Time via Topology 


We have seen in the last section that if one considers all the topological information encoded in the microstruc¬ 
ture evolution of a solution path of the Cahn-Hilliard-Cook model then one can use the classification 
schemes (CA) and (Ck) to recover the total mass ju with high accuracy. Of course the total mass is only 
one of several parameters in (Q which affect the topology, and it would be interesting to also determine how 
strong the influence of these other parameters on the topology is. 

Rather than pursuing this line of inquiry, we now consider a different question. For this, consider again the 
simulations of the Cahn-Hilliard-Cook model, but now for a fixed and a-priori known mass value ju. Then a 

simulation of @ produces solution snapshots along M equidistant timesteps fi. tM, and for each time ti, 

computational topology can be used to create the associated persistence landscapes in dimensions 0 and 1. 
Suppose now that we pick up one of these persistence landscapes at random. Can we recover the actual 
time ti which corresponds to this landscape purely from the encoded topological information? In other words, 
for fixed mass value n, is it possible to determine the time at which a solution snapshot has been taken, 
purely based on topological information? We will see in the following that this can in fact be done, and in 
many cases with surprising accuracy. 

As we saw in the previous section, our classification technique can only recover mass information up to 
a certain tolerance, and the situation is no different in the snapshot time case. Therefore, we only use the 
data set (D3) in the present section, which produces M = 25 times t^ in the interval (0,Te]. Recall that 

this data set considers the N = 51 mass values p, = 0.01 /? for n = 0.50, and for each mass value one 

has R = 100 stochastically independent repetitions of the simulation. 

Based on the composition of the data set (D3), we apply the classification schemes (CO) and (Cl) from 

Section [rll for every fixed mass value p, to the snapshot time classes Ci. Cm, which correspond to the M 

snapshot times ti . tM in increasing order. Then we proceed as follows. 


Based on Rj < R = 100 training runs each from the M = 25 considered snapshot times ti, we 
determine the averaged classifier of the class in dimension k as defined in (14) with S = Rt 
and n = m = 1 . M. 


• We then try to classify each of the remaining {R — Rt)M simulations in the data set for fixed mass 
value p, by either using the classification scheme (CO) or (Cl). 


As before, in each case we note the index m obtained from the classification scheme, and compare it 
to the snapshot time index mtrue of the underlying simulation. If m = mtrue, then we have obtained a 
“Hit”, otherwise the classification “missed by \m — mtrue\'- Clearly the integer value \m — mtrue\ is directly 
proportional to the time difference incurred by the classification, with the proportionality factor being given 
by Te/M = 0.04 Te. 
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Figure 11: Recovering the snapshot time from persistence information using the L^-norm and training set 
size Rt = 50. Each panel contains a heat map, whose horizontal axis corresponds to the A/ = 51 mass 
values n, and the vertical axis to the M = 25 snapshot times, in increasing order. Color indicates frequency 
of observation, as indicated by the colorbar. Panels in the left column are for scheme (CO), the right column 
corresponds to (Cl). From top to bottom, the three rows show the likelihoods of exact hits, misses by 
exactly one, and misses by at least two. 


Rather than performing an exhaustive study as in the previous section, we use the already obtained insight 
to focus our approach. Since the /.°°-norm classification consistently performed significantly worse than the 
L^- and -norm classifications, we only run simulations for the latter two. Flowever, in contrast to the 
experiments presented in Section 3^ it turns out that for recovering the snapshot time both (CO) and (Cl) 
are superior to (CA). For this reason, we only consider the dimension-dependent classification schemes below. 

For the case of L^-norm classifications using the schemes (CO) and (Cl), and for training set size Rt = 50, 
the results can be found in Figure [1^ For the sake of brevity, we have presented the data in the form of 
heat maps, rather than providing precise tables. Each panel in the figure represents frequencies which are 
indicated by colors, and which correspond to the ranges of the adjacent colorbars. The horizontal axes 
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Figure 12: Recovering the snapshot time from persistence information using the L^-norm and training set 
size Rt = 50. Each panel contains a heat map, whose horizontal axis corresponds to the A/ = 51 mass 
values n, and the vertical axis to the M = 25 snapshot times, in increasing order. Color indicates frequency 
of observation, as indicated by the colorbar. Panels in the left column are for scheme (CO), the right column 
corresponds to (Cl). From top to bottom, the three rows show the likelihoods of exact hits, misses by 
exactly one, and misses by at least two. 


correspond to the N = 51 mass values while the vertical axes are for the M = 25 snapshot times, both 
in increasing order. Panels in the left column are for classification scheme (CO), while the right column 
corresponds to (Cl). From top to bottom, the three rows show the likelihoods of exact hits, misses by 
exactiy one, and misses by at ieast two. Analogous results for the case of the /.^-norm are depicted in 
Figure[I^ In each of the cases we have cross-validated the results by choosing different training sets of the 
same size, and the results are quantitatively similar. Notice that in this setting, for each mass value n we 
are classifying (/? — Rt)M = 1250 microstructures according to their snapshot time class. 


A first glance at the heat maps in Figures and especially the ones in the respective bottom rows, 
confirms our statement from the beginning of this section. Averaged persistence landscapes can indeed be 
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used to accurately deduce the snapshot time at which a given microstructure occurs. In most cases, the 
match is an exact hit, while in a few cases one is off by one index. In more detail, one can draw the following 
conclusions from the data in these figures: 

• For every considered mass value fj,, there exists an interval of the form [a^,l3^] c (0, Tg] over which 
the snapshot time can be determined exactly with probability close to one. 

• The times > 0 are small, and more or less constant with Only at the beginning of the simulations 
does it seem to be difficult to pin down the exact snapshot time. The volatile timeframe (0, a^) corre¬ 
sponds to the initial rapid smoothing regime in the parabolic stochastic partial differential equation Q, 
during with the high-frequency terms in the random initial conditions decay due to dissipativity. 

• The value of/3^ is close to about 60% of the endtime Tg for mass values ju close to zero, and decreases 
as jLi increases. At the final mass value iJ, = 0.5 one has Pfj, « 0.3 Tg. For larger simulation times, i.e., 
times in the interval Tg], it appears to be more difficult to precisely locate the snapshot time, and 
this becomes more pronounced as the mass n increases. This timeframe is well into the coarsening 
regime, where microstructures, particularly of droplet type, only move slightly, and no longer change 
their topology frequently. Thus, the classification results are not surprising. 

Notice that in all panels of the two figures, the probabilities seem to jump significantly along the topmost 
line. In other words, the classification results for the time tM are significantly different from the one for t/w-i. 
This is due to the fact that for all interior times ti < Tg one can miss by one by either classifying the time 
as ti-i or while at the last time, anything that would have been misclassified as a larger time will most 
likely be classified as tM- 

For the above results we used different norms, but kept the training set size fixed at Rj = 50. In order to 
see how this size affects the results, we have repeated the simulations, but this time with smaller training sets 
of size Rt = 20. The corresponding results are shown in Figures [T^and M Notice that despite the smaller 
training set size, and the thereby induced larger simulation size of {R— Rt)M = 2000 microstructures which 
have to be classified, the results are almost the same. This indicates again that the topological information 
which is captured by the random microstructures and their evolution carries enormous information, enough to 
accurately detect the snapshot time, and therefore the stage of material decomposition. The only apparent 
exception to this is the case of large mass values ju « 0.5 and time values well into the coarsening regime. In 
the first case, the microstructures are of moving droplet type as shown in the bottom row of Figure [TJ and 
in the latter case microstructure evolve slowly since there is less energy to dissipate, and therefore there are 
fewer topology changes. But during the initial spinodal decomposition regime, and well into the coarsening 
regime, topological information can be used to find the snapshot time. 

4 Summary and Conclusions 

Homology is a way of quantifying complicated topological spaces through a sequence of discrete objects, in 
fact, through a sequence of integers in its most reduced form. This information is invariant under continuous 
transformations of the underlying space, and in particular, it does not contain any size or “metric" information 
of the considered object. The only structures that "count", in the true sense of the word, are /c-dimensional 
holes, which cannot be contracted within the space. In this form, homology has long been used in classical 
mathematics to determine when two spaces are different within a certain category, be it that they are not 
homeomorphic or not homotopy equivalent. For data analysis, however, often the opposite question is of 
interest. Can homology be used to say that two objects are similar, and if so, how similar are they? 

At its core, the above question is concerned with the amount of information that is retained by homology. 
In the present paper, we study this problem in a well-defined and important situation, namely phase separation 
dynamics in binary metal alloys as described by the stochastic Cahn-Hilliard-Cook model Q. The equation 
can be used to create realistic evolutions of phase separating microstructures, which due to their inherent 
complexity cannot easily be classified using standard techniques. By using persistence landscapes, we could 


28 



Figure 13: Recovering the snapshot time from persistence information using the L^-norm and training set 
size Rt = 20. Each panel contains a heat map, whose horizontal axis corresponds to the A/ = 51 mass 
values n, and the vertical axis to the M = 25 snapshot times, in increasing order. Color indicates frequency 
of observation, as indicated by the colorbar. Panels in the left column are for scheme (CO), the right column 
corresponds to (Cl). From top to bottom, the three rows show the likelihoods of exact hits, misses by 
exactly one, and misses by at least two. 

obtain a sequence of discrete objects, which encapsulates the topology evolution of a solution path. Since 
the Cahn-Hilliard-Cook model is stochastic in nature, and in addition, since the initial alloy state is never 
known precisely, one can only hope to capture typical behavior, and we do this by averaging these persistence 
landscapes. It was demonstrated in this paper, that these averaged landscapes encode enough information 
to make the following decision problems solvable: 

• Given only topological information about a specific microstructure evolution, can we determine the 
total mass n of the underlying experiment? 

• Given only topological information about a specific microstructure at a given mass value n, can we 
determine the snapshot time during the decomposition process at which this pattern was observed? 
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Figure 14: Recovering the snapshot time from persistence information using the L^-norm and training set 
size Rt = 20. Each panel contains a heat map, whose horizontal axis corresponds to the A/ = 51 mass 
values n, and the vertical axis to the M = 25 snapshot times, in increasing order. Color indicates frequency 
of observation, as indicated by the colorbar. Panels in the left column are for scheme (CO), the right column 
corresponds to (Cl). From top to bottom, the three rows show the likelihoods of exact hits, misses by 
exactly one, and misses by at least two. 

These questions can be answered by comparing the provided topological information to the nearest averaged 
classifiers, which are created through training sets. 

We would like to stress that we did not set out to develop a method which determines the total mass ju 
from topological information — that can obviously be done easier by just finding the average gray level of 
any solution snapshot. Rather, we wanted to demonstrate that the seemingly stark reduction in information, 
which happens during the passage from a complicated microstructure to its persistent homology, still retains 
an incredible amount of information about the underlying dynamical process. Yet, due to its inherently 
discrete nature, this topological approach shows considerable promise for analyzing and relating experimental 
data to numerical models, and this avenue has to be explored further in future work. 
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